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EFFICIENT  DISCRETE  FOURIER  REPRESENTATION 
OF  PULSE  RESPONSES 

INTRODUCTION 


The  problem  considered  in  this  report  is  the  modeling  of  the  transient  response 
of  a  linear  system  to  an  incident  pulse  (finite  time  series).  In  mathematical  modeling, 
this  problem  is  solved1  through  the  convolution  of  a  time  domain  pulse  with  a  model¬ 
generated  transfer  function  (spectral  Green’s  function)  of  the  linear  system  under 
consideration.  This  report  is  concerned  with  a  particular  class  of  linear  systems,  i.e.,  those 
that  are  separable  into  a  superposition  of  multiple  linear  subsystems.  Such  problems  arise 
often  in  high-frequency  acoustics  applications,  where  ray  theory  is  applied  and  the  linear 
acoustic  system  is  represented  as  a  superposition  of  many  acoustic  rays  of  energy. 

The  convolution  of  a  time  domain  pulse  with  a  linear  system  is  a  well-known 
problem  in  linear  systems  theory2,3  that  arises  often  in  areas  such  as  signal  processing  and 
control  systems.  The  techniques  used  are  based  on  a  Fourier  representation  of  the  system 
and  the  incident  pulse  in  the  frequency  domain.  Simple  convolution  provides  the  system 
output  in  the  time  domain.  The  superposition  principle  allows  addition  of  subsystem 
responses  in  either  the  time  or  frequency  domains,  the  latter  of  which  is  much  more 
straightforward  yet  computationally  inefficient.  By  examining  the  practical  nuances  of 
the  time  domain  approach,  a  computationally  efficient  method  of  superposition  of 
subsystems  can  be  developed. 

In  this  report,  a  numerically  efficient  technique  is  presented  for  the  computational 
problem  of  mathematically  modeling  the  convolution  of  a  transient  pulse  with  a  set  of 
subsystem  transfer  functions.  Both  the  theoretical  and  practical  limitations  of  this 
technique,  as  well  as  the  standard  approach,  are  presented.  The  results  of  the  technique 
are  shown  in  the  context  of  two  problems  of  Navy  importance:  acoustic  scattering  from  a 
multiple-highlight  structure  and  acoustic  propagation  modeling  using  eigenray  bundles. 
Results  are  presented  in  the  context  of  both  computational  efficiency  as  well  as  numerical 
accuracy.  The  technique  is  developed  in  a  general  framework  such  that  its  extension  to 
other  applications  is  straightforward.  The  report  concludes  with  general  guidance  on  the 
utility  of  the  technique  in  applications. 
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TRANSIENT  PULSE  RESPONSE  MODELING 
FOR  A  SINGLE  SYSTEM 


The  modeling  of  steady-state  phenomena  is  a  common  practice  in  many  areas  of 
science  and  engineering.  Often,  these  are  the  only  models  that  can  be  analytically  and/or 
numerically  solved  in  an  efficient  manner;  therefore,  many  engineering  decisions  are 
made  based  on  steady-state  models  even  when  the  process  under  consideration  is  not 
steady-state.  However,  in  the  case  of  linear  input/output  systems,  the  steady-state 
response  can  be  used  to  determine  transient  behavior.  In  this  section,  the  standard 
approach  for  modeling  transient  responses  using  steady-state  models  is  developed  and 
discussed. 

The  process  that  is  being  modeled  is  assumed  to  have  a  single  (scalar)  input  and 
single  (scalar)  output.  This  is  an  assumption  that  is  made  for  notational  convenience  only, 
as  the  multi-input,  multi-output  case  is  a  simple  extension  with  no  theoretical  differences. 
Furthermore,  assume  that  the  process  being  modeled  is  representable  by  an  impulse 
response  function  (or  Green’s  function)  given  by  h(t).  Using  the  theory  of  linear 
systems,1,2  the  response  y(t )  of  the  system  to  a  pulse  u(t)  is  given  by  the  convolution 


y(t)  =  f  h(t)u(t  -t)dr. 

J—oo 


(1) 


The  only  assumptions  required  in  equation  (1)  are  that  the  system  is  a  linear  time- 
invariant  system  and  the  input  pulse  u(t )  is  a  real  excitation.  Note  that  if  the  system  is 
causal,  then  h(t)  =  0  for  all  t  <  0;  thus,  in  that  case  the  integral  may  be  taken  from  0  to 
+°° .  The  expression  in  equation  (1)  can  be  placed  in  the  frequency  domain  by  Fourier 
transforming  both  sides  to  yield 

Y(co)  =  H(cq)U(cq),  (2) 

which  is  a  standard  result  from  elementary  linear  systems  theory.  In  expression  (2),  the 
term  Y(cd)  is  used  to  represent  the  Fourier  transform  of  y(t),  and  similarly  for  H{oS)  and 
U(co).  In  this  notation,  H(cd)  is  referred  to  as  the  transfer  function  of  the  system. 

If  the  system  under  examination  is  subjected  to  a  unity  amplitude  steady-state  input 
of  frequency  (Do ,  then  the  resulting  output  is  given  by  (upon  inspection  of  equation  (2)) 

Y(CQ0)  =  H(cq0).  (3) 

Thus,  the  transfer  function  H(o)o)  at  any  frequency  coq  can  be  found  as  the  response  of  the 
system  to  a  steady-state  input  of  frequency  ay>.  This  provides  a  convenient  means  of 
numerically  determining  a  transfer  function  by  running  a  steady-state  model  at  any 
required  frequency.  Since  many  mathematical  models  of  complex  systems  are  derived  in 
a  steady-state  manner,  this  is  usually  a  straightforward  extension  of  an  existing  model.  It 
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should  be  noted,  however,  that  this  technique  only  provides  samples  of  the  transfer 
function  at  a  discrete  set  of  frequencies,  rather  than  the  continuous  function  H{oS). 
Furthermore,  the  transfer  function  is,  in  general,  a  complex  quantity,  and  thus  any  model 
must  accurately  produce  both  the  magnitude  and  phase  of  the  steady-state  output  for 
application  as  a  sampled  transfer  function.  Many  existing  steady-state  models  are  only 
validated  for  magnitude,  and  phase  validation  over  the  frequency  range  of  application  is  a 
necessity  for  transient  response  modeling.  Thus,  additional  validation  of  the  steady-state 
model  is  often  required  prior  to  its  use  in  transient  response  modeling. 

The  mathematical  modeling  problem  of  interest  is  to  convolve  a  given  (non- 
steady-state)  pulse  with  a  system  defined  by  an  impulse  response  function.  For  the 
purposes  of  this  report,  it  will  be  assumed  that  the  system  in  question  has  been  adequately 
modeled  for  steady-state  response  and  that  a  numerical  model  of  that  process  exists.  It  is 
further  assumed  that  the  steady-state  model  has  been  validated  for  both  magnitude  and 
phase  over  the  frequency  range  of  interest.  Some  examples  of  available  models  in 
acoustics  are  finite-element  models  of  structural  acoustic  response,  environmental 
propagation  models  for  a  fixed  source/receiver  pair,  and  acoustic  scattering  models  for  a 
single  highlight. 

Let  h(t)  denote  the  system’s  impulse  response  function  (Green’s  function). 
Accordingly,  let  H(O)j)  denote  the  transfer  function  of  the  system  (steady-state  model 
output)  at  a  frequency  C0j,  which  is  a  sampling  of  the  Fourier  transform  of  h(f)  at  a 
specific  frequency.  Let  p(t)  represent  a  real  input  pulse.  It  is  assumed  that  p(t)  is  non¬ 
zero  only  over  a  finite  interval  of  time  from  t  =  to  to  t  =  t\.  From  expressions  (1)  and  (2), 
it  is  easily  seen  that  the  system  response  is  given  by 

y(0  =  F"1{r(G))} 

=  F-1{tf(<o)P(<y)}  (4) 

=  F-1{H(o))F{p(0}}, 


where  F{ }  represents  the  Fourier  transform  and  F_1{ }  represents  the  inverse  Fourier 
transform.  Using  a  Fourier  transform  pair  of  the  form 


G{.(0)=rg{t)e~imdt, 

•f— oo 

(5a) 

and 

g{t)  =  -^-V~G((Q)eimd(D, 

In  J-~ 

(5b) 

it  can  be  shown  that  equation  (4)  is  equivalent  to 

y(t)  =  ^  j[J  H(co)p(t)em'-r)dTdco. 

(6) 
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Note  that  any  valid  Fourier  transform  pair  will  produce  the  same  result  as  in  equation  (6). 
The  expression  in  equation  (6)  will  be  referred  to  as  the  general  form  of  the  system 
response  function  for  a  given  input.  Because  p(t)  is  defined  to  be  nonzero  only  over  a 
finite  interval  of  time,  equation  (6)  can  be  rewritten  in  the  form 


y(0  = 


(  'L 


*  J-oo 


dco, 


(7) 


which  makes  explicit  the  finite  time  duration  of  the  input  pulse. 

The  form  of  equation  (7)  is  a  theoretical  result  on  the  use  of  steady-state  system 
behavior  to  determine  a  response  due  to  an  arbitrary  finite  time  duration  pulse.  In 
practice,  the  values  on  the  right  side  of  this  equation  are  not  available  in  a  continuous 
form;  they  are  only  sampled  values.  Furthermore,  in  practice,  only  the  real  part  of  y(t)  is 
needed,  and,  in  the  remainder  of  this  report,  it  is  assumed  that  the  system  response  is 
given  by  the  real  part  of  y(t).  When  the  transfer  function  is  sampled  at  a  discrete  set  of 
frequencies  {cty},  the  frequency  integral  can  be  replaced  by  a  summation  to  obtain 


y(r)« 


2  7t  ) 


^Hico^dcoj  p(r)eiWj('~T)dT 


(8) 


where  5&>is  the  uniform  spacing  between  frequency  samples  {to,}.  By  further  assuming 
that  the  input  pulse  is  only  available  in  the  form  of  a  sampled  time  series  sampled  at  the 
discrete  set  of  time  points  {4},  the  system  response  function  becomes 


y(t)  ~ 


1 


2n{tx  -4) 


Re 


)8co^  p(tk  )e,<0j('  Ti)5t 


(9) 


where  <5ris  the  uniform  spacing  between  time  samples.  This  expression  is  a  discrete  form 
of  equation  (7)  that  makes  use  of  the  finite  time  extent  of  the  incident  pulse  p(t).  The 
necessary  values  of  these  spacings  for  an  accurate  solution,  as  well  as  the  range  of 
frequency  samples  required,  are  issues  that  are  addressed  using  sampling  theory,  along 
with  some  basic  knowledge  of  the  system  being  modeled.  It  should  be  noted  that  fast 
Fourier  transform  techniques4  can  be  used  to  represent  the  integral  found  in  equation  (8). 
Such  techniques  will  improve  overall  computational  efficiency,  but  do  not  affect  the 
results  of  the  comparative  study  of  two  techniques  that  are  examined  in  this  report. 

The  sampling  theorem  can  be  summarized  as  follows.  A  time  signal  whose 
spectrum  (Fourier  representation)  has  only  zero  components  outside  the  finite  frequency 
band  ffom/=/i  to/ =/2  can  be  uniquely  determined  by  its  values  at  (uniform)  time 
samples  spaced  by 


4 


(10) 


A t- 


1 

2C/2-/,)’ 


The  sampling  interval  given  by  equation  (10)  is  referred  to  as  the  Nyquist  sampling  rate. 
This  is  only  a  theoretical  result,  since  all  real  pulses  (i.e.,  those  with  finite  time  extent) 
must  have  frequency  components  over  an  infinite  frequency  band.  However,  practice 
shows  that,  for  acoustics  applications,  a  band  that  encompasses  those  frequency 
components  whose  Fourier  magnitudes  are  greater  than  one-tenth  of  the  maximum  spectral 
magnitude  provides  an  acceptable  frequency  band.  This  limitation  on  the  practical 
frequency  domain  of  the  input  pulse  provides  the  necessary  time  sampling  of  the  pulse  (via 
equation  (10))  as  well  as  the  time  sampling  of  the  resulting  output  time  series,  which  are 
sampled  at  the  same  rate.  Another  practical  rule-of-thumb  in  the  use  of  sampling  to  create 
digital  pulses  is  to  sample  at  twice  the  Nyquist  rate,  or  at  a  At  that  is  one-half  of  what  is 
given  by  equation  (10). 


A  second  issue  that  arises  in  the  direct  use  of  the  sampling  theorem  is  that  the 
system  must  be  basebanded  to  directly  apply  the  sampling  theorem  to  a  frequency  range 
that  does  not  extend  to  zero.  While  basebanding  may  be  a  common  practice  in  signal 
processing,  it  is  difficult  to  perform  in  modeling  since  both  the  input  pulse  and  the  system 
impulse  response  must  be  basebanded  in  an  identical  manner  in  order  to  arrive  at  a  result 
that  can  be  transformed  back  to  the  proper  frequency  range.  To  avoid  this  problem,  the 
system  will  be  sampled  assuming  that  the  lower  limit  of  the  frequency  range  (ft  in  equation 
(10))  is  equal  to  zero.  This,  along  with  the  other  practical  sampling  issues  stated  above, 
leads  to  a  practical  time  sampling  rule  of 


At  = 


4/2 


where 


/2  =max{/}3  \ 


max 


Fjpi  t) 


10 


)  J 


(Ha) 

(Hb) 


Expression  (11)  may  look  complicated  in  mathematical  notation,  but  it  merely  states  that 
the  sampling  frequency  is  four  times  the  maximum  frequency  of  a  frequency  band  given 
by  the  20-dB  down  points  of  the  power  spectrum  of  the  incident  pulse. 

The  frequency  limitations  used  in  the  development  of  the  time  sampling  in 
equation  (1 1)  can  also  be  used  to  develop  a  compatible  limitation  on  the  frequency  band 
bounds  for  the  transfer  function  H(co).  Define /1  by 


fi  =  min(/)  3 


|F{pfr)}|  < 


^max|F{p(Qp 

10 


(12) 
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Now/i  and/2  define  the  bounds  for  the  frequency  band  for  the  transfer  function  evalua¬ 
tion  of  the  system  response.  The  remaining  parameter  necessary  for  the  evaluation  of  the 
system  response  function  given  in  equation  (9)  is  to  determine  the  frequency  step  size  to 
use  in  establishing  the  set  {C0j}. 

The  step  size  to  use  in  frequency  stepping  is  an  issue  that  is  different  for  modeling 
purposes  than  it  is  in  the  signal  processing  community.  The  frequency  step  that  is  used 
affects  the  inverse  Fourier  transform  (equation  (5b))  by  forcing  a  periodicity  of  the 
resulting  time  series.  That  is,  if  an  inverse  Fourier  transform  is  made  in  a  discrete  manner 
with  steps  of  8(0,  then  the  continuous  inverse  Fourier  transform  becomes  the  discrete 
inverse  Fourier  transform  as  follows: 

1  400 

*«)  =  =>  g(j)=  ■  (13) 


From  the  above  expression,  it  is  easily  seen  that  any  time  shift  of  the  form  t^t  +  T, 
where  T  is  given  by 


T 


f  2n\ 


(14) 


for  any  integer  m,  will  not  change  the  value  of  g(t).  Thus,  any  discrete  inverse  Fourier 
transform  that  is  sampled  at  a  frequency  step  of  Sco  will  have,  by  necessity,  a  consequence 
of  periodicity  in  time  with  a  period  of  T  as  given  in  equation  (14).  For  the  application  of 
mathematical  modeling  of  linear  input/output  systems  as  addressed  in  this  report,  this 
implies  that  the  model  result  y(t)  will  be  periodic  with  a  period  T  as  given  in  equation  (14). 
Such  a  periodicity  places  a  constraint  on  the  maximal  value  to  use  for  Sco,  which  is  simply 

-  .  lit 

5co<—,  (15) 

1o 


where  To  represents  the  minimum  length  of  time  for  which  there  is  output  from  the  input/ 
output  process.  This  time  extent  must  be  obtained  from  the  physics  of  the  process  being 
modeled,  and  is  not  purely  a  mathematically  derived  quantity  for  a  process  modeling 
problem.  The  practical  limitation  on  this  value  for  acoustic  applications  is  what  makes  the 
direct  method  of  Fourier  realization  impractical  for  many  acoustic  modeling  tasks. 

The  number  of  points  needed  in  time  and  frequency  to  properly  construct  an  output 
time  series  can  be  quantified  using  the  results  to  this  point.  Using  equations  (11)  and  (12), 
along  with  the  input  pulse,  one  can  determine  the  frequency  range  and  time  step  needed. 
Furthermore,  the  pulse  length  itself  gives  the  extent  of  time  needed  to  be  modeled  for  only 
the  input;  the  time  extent  of  the  output  will,  in  general,  be  longer  than  that  of  the  input  for 
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transient  system  responses.  Based  on  the  output  time  extent,  equation  (15)  can  be  used  to 
set  a  frequency  step.  Thus,  the  number  of  frequency  points  m  required  is  given  by 


a  RA-ZlRI.  (16) 

where  /  and/2  are  as  in  equations  (11)  and  (12).  The  notation  [x~|  is  a  common  practice 
in  mathematics  literature  used  to  represent  the  smallest  integer  greater  than  or  equal  to  x. 
Similarly,  the  symbol  j_xj  represents  the  largest  integer  less  than  or  equal  to  x.  The 
number  of  time  points  required  for  the  input  nin  and  the  output  nout  are  given  by 

-FVa-fc-Ol.  (17) 


= 


h _ h 

At 


and 


n 


out 


(18) 


Therefore,  the  double  sum  in  equation  (9)  must  be  taken  at  m  frequency  points 
and  riin  time  points  for  a  total  of  M  time-frequency  points,  where 


M  -m-nin=  [4  T0f2(  f2  -  /  )( tx  - 10  )]. 


(19) 


Each  of  these  M  time-frequency  points  are  used  in  the  evaluation  of  output  at 
each  of  the  nou,  output  time  steps.  From  expression  (19),  it  is  obvious  that  M  is  large 
when  the  input  signal  has  a  large  time-bandwidth  product  and/or  when  the  output  time 
extent  required  is  very  long.  When  both  of  these  conditions  occur,  the  problem  becomes 
computationally  very  costly. 
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TRANSIENT  PULSE  RESPONSE  MODELING 
FOR  MULTIPLE  SUBSYSTEMS 


When  a  system  is  composed  of  multiple  independent  subsystems,  the  practical 
application  of  transient  pulse  response  modeling  becomes  computationally  costly.  This 
problem  arises  frequently  in  acoustic  modeling,  such  as  propagation  modeling  using  a 
summation  of  acoustic  eigenrays,  or  acoustic  scattering  from  a  structure  composed  of 
multiple  independent  highlights.  In  these  cases,  steady-state  models  of  the  individual 
subsystems  have  been  created  because  of  computational  and/or  theoretical  simplicity.  Due 
to  the  linearity  of  acoustics,  the  results  of  these  subsystems  may  be  combined  using  linear 
superposition  to  arrive  at  a  single  steady-state  result  for  the  composite  system.  In  this 
section,  the  extension  of  the  single  system  case  to  multiple  subsystems  is  explained,  and  an 
efficient  computational  method  for  handling  the  multiple  subsystem  case  is  derived. 

Assume  that  a  system  with  transfer  function  H(co)  is  composed  of  N  independent 
subsystems,  with  linear  transfer  functions  {#„(&>)}  such  that 

ff  (a) -£»„(©).  (20) 

«=1 

This  kind  of  separation  into  subsystems  is  often  done  when  steady-state  models  (analytical 
and/or  computational)  of  the  individual  subsystems  are  readily  available.  However,  in 
order  to  apply  equation  (20)  to  a  physical  modeling  problem,  a  number  of  conditions  must 
be  met.  First,  each  of  the  subsystems  must  be  a  linear,  time-invariant  system.  The  time- 
invariance  implies  that  an  input  of  u{t  +  td),  for  any  time  delay  td,  will  produce  the  output 
y(t  +  td).  Note  that  this  type  of  time-invariance  will  maintain  any  relative  time  delays 
between  different  subsystems,  since  the  delay  is  applied  to  the  input,  not  the  system 
impulse  response  function.  Secondly,  the  subsystems  must  be  linearly  independent  from 
one  another  to  allow  the  superposition  principle  to  hold.  Finally,  each  of  the  subsystems 
must  be  causal,  i.e.,  there  can  be  no  dependence  on  any  past  input  or  state  for  any  of  the 
given  subsystems.  Thus,  the  response  to  any  of  the  given  subsystems  does  not  depend  on 
any  of  the  other  subsystems. 

The  standard  approach  to  modeling  a  multiple  subsystems  problem  is  to  use 
equation  (20)  directly  by  summing  the  component  transfer  functions  in  the  frequency 
domain  and  using  that  result  as  the  composite  system  transfer  function.  This  result  is 
accurate,  but  may  be  extremely  computationally  costly.  Consider  the  following  simple 
example.  A  system  of  N  subsystems  has  transfer  functions  that  are  given  by 


ff„(o»  =  V*. 


(21) 


for  positive  A„  and  The  time-domain  output  of  the  composite  system  is  given  by 
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(22) 


y(t)  =  '£Anu(t+%n), 

n= 1 

which  is  a  sequence  of  scaled  and  time-shifted  replicas  of  the  input  pulse.  If  the  input 
pulse  has  a  time  extent  that  is  nonzero  over  a  range  from  to  to  t\,  then  the  resulting 
necessary  time  extent  of  the  output  will  be  given  by 

T0  =  (max{£„}- min{£n  })+(f,  -t0).  (23) 

This  number  directly  affects  the  number  of  frequency  samples  needed  for  accurate  pulse 
output  modeling  as  shown  by  equation  (16).  The  form  of  the  expression  in  equation  (23) 
shows  that  the  number  of  frequency  samples  required  can  (in  this  particular  case)  scale 
linearly  with  the  difference  in  time  delays  between  subsystems. 

The  result  of  this  simple  example  shows  that  the  number  of  frequency  samples 
required  can  become  very  large  when  the  relative  time  delays  between  subsystems 
become  large.  In  this  example,  the  number  of  frequency  samples  required  for  the  output 
representation  is  found  from  equation  (16)  to  be 

m  =  (f 2  -  fi  ){h -to)+(f2  -  fi  )(max{|n }  -  min{|„ }).  (24) 

The  first  term  in  this  expression  is  equal  to  the  number  of  frequencies  required  to  model 
the  input  pulse  in  a  Fourier  representation,  which  is  added  to  a  term  that  accounts  for  the 
different  time  delays  found  in  the  various  subsystems.  Most  signal  processing  applica¬ 
tions  only  window  a  time  series  over  the  length  of  the  input  pulse;  therefore,  in  that  case, 
the  second  term  in  this  expression  is  negligible.  However,  in  the  mathematical  modeling 
of  pulse  responses,  both  terms  in  equation  (24)  are  required  for  a  satisfactory  solution. 
This  use  of  the  direct  summation  of  subsystem  transfer  functions  for  modeling  will  be 
referred  to  as  the  direct  approach  for  the  remainder  of  this  report. 

In  the  general  setting,  consider  a  system  comprised  of  N  subsystems  with  transfer 
functions  {Hn  (ft))}.  Each  of  the  subsystem  responses  can  be  placed  in  the  form  of  a 
secondary  transfer  function  Hn  (tu)  and  a  time  delay  of  such  that  the  subsystem 
transfer  function  is 


Hn(co)  =  Hn(co)e 


ic°Zn 


(25) 


where  the  secondary  transfer  function  contains  no  additional  time  delays.  That  is,  if  the 
input  pulse  begins  at  time  t  =  0,  then  the  secondary  transfer  function  has  an  output  that 
also  begins  at  time  ?  =  0.  In  this  way,  the  time  delay  feature  of  the  subsystem  transfer 
function  has  been  separated  from  the  remaining  effects  of  the  subsystem  transfer 
function. 


9 


Once  all  of  the  subsystem  transfer  functions  have  been  properly  separated  as  shown 
in  equation  (25),  the  result  of  the  total  system  output  can  be  derived  from  equation  (9)  as 


y(0» 


1 


27r(tj  t0  ) 


N 

2  Re 


n~\ 


X  H,  (ffl ,  )e“'M-  ’&>£  P(t*  y1’*  Sr 


(26) 


This  expression  is  even  more  computationally  intensive  than  the  transfer  function 
summation  method  (the  direct  approach);  however,  equation  (26)  is  greatly  simplified 
under  a  constraint.  If  the  subsystem  time-delays  are  all  constrained  to  lie  on  steps  of 
the  output  time  axis  t,  then  the  exponential  expression  in  the  middle  of  equation  (26)  is 
just  a  discrete  time  shift  of  the  axis  of  the  form 

**-**  +  £,.  (27) 

for  each  of  the  N  subsystems.  Note  that  if  the  added  constraint  were  not  made,  the 
evaluation  of  equation  (26)  along  a  discrete  time  axis  would  be  impractical,  since  each 
subsystem  would,  in  general,  be  on  a  slightly  shifted  time  axis.  With  the  added  constraint, 
the  final  form  can  be  written  as 

y(r) «  £  Re{y„  (*  +  £,)},  (28) 


where 


M)  = 


2  it 


(29) 


and  P{(oj)  is  the  discrete  form  of  the  Fourier  transform  (as  in  equation  (5))  of  pit)  evaluated 
at  a  frequency  of  CDj.  The  use  of  equations  (25),  (28),  and  (29)  combined  creates  a  computa¬ 
tionally  efficient  method  of  representing  the  combination  of  multiple  subsystem  responses. 

In  practice,  the  assumption  that  all  of  the  subsystem  delays  (as  shown  in  equation 
(25))  lie  exactly  on  the  discrete  steps  of  the  output  time  axis  is  impractical.  In  that  case, 
the  subsystem  transfer  function  is  first  placed  in  a  form  similar  to  equation  (25),  which  is 
given  by 

=  (30) 

where  the  total  time  delay  of  the  subsystem  is  A tn  =  (£n  +  £„).  Since  the  model  output  for 
each  subsystem  will  contain  the  total  time  delay,  the  parameters  and  £„  are  somewhat 
arbitrary.  Any  combination  whose  sum  equals  the  total  subsystem  time  delay  will  be 
accurate,  but  computational  efficiency  is  met  when  they  are  chosen  such  that 
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At. 


and  Cl=Arn-|B, 


(31) 


where  tstep  is  the  uniform  spacing  of  samples  in  the  output  time  axis.  Now  the  transfer 
function  representation  for  the  subsystem  is  given  as  shown  in  equation  (25)  with 

H,(a)  =  H„(<oV“-.  (32) 

In  this  manner,  the  transfer  functions  of  each  of  the  subsystems  have  been  accurately 
shifted  so  that  the  subsystem  time  delays  all  line  up  with  the  time  steps  on  the  output  time 
axis.  By  performing  this  transformation  for  the  slight  shifting  along  the  time  axis,  the 
interference  effects  (both  constructive  and  destructive)  that  occur  due  to  multiple 
subsystems  overlapping  in  time  are  maintained  exactly.  The  outputs  from  the  subsystems 
can  then  be  added  together  coherently  on  the  output  time  axis  to  obtain  the  same  result 
(up  to  round-off  error)  that  would  be  obtained  using  a  single  transfer  function  as  in 
equation  (20).  The  use  of  equation  (31)  eliminates  any  errors  due  to  slight  shifting  of  the 
subsystems  along  the  time  axis  to  line  up  with  the  time  steps. 

The  computational  requirements  of  a  system  modeled  using  the  efficient  technique 
can  be  addressed  in  the  same  manner  as  for  the  direct  approach.  Let  T„  denote  the  time 
extent  of  each  of  the  output  pulses  for  the  subsystems.  In  the  case  of  an  input  pulse  that 
exists  from  t  =  t0  to  t  =  t\,  the  output  pulse  for  the  n*  subsystem  exists  (is  nonzero)  over 
the  time  interval  from  t  =  (t0  +  ^n)  to  t  =  (t0  +  %„)  +  T„.  From  equation  (29),  it  is  seen  that 
the  output  pulses  for  each  subsystem  will  be  computed  independently;  thus  all  of  the  Tn 
values  are  independent.  Each  of  the  subsystems  has  its  own  frequency  summation  in 
equation  (29);  therefore,  using  equation  (16),  the  number  of  frequency  samples  that  must 
be  computed  for  each  subsystem  is 

(33) 

It  must  be  remembered  that  the  number  of  frequency  samples  required  is  also  the  number 
of  times  that  the  steady-state  model  must  be  exercised  to  obtain  the  transfer  function 
samples  {Hn(C0j)}.  In  the  case  of  the  efficient  technique,  each  subsystem  only  is  run  at  mn 
frequencies  (as  shown  in  equation  (33)),  while  the  direct  approach  requires  all  subsystems 
to  be  computed  at  the  m  frequency  samples  as  given  by  equation  (16).  Assuming  that  all 
subsystems  require  the  same  amount  of  computational  effort  per  frequency,  and 
representing  the  effort  for  computing  the  transfer  function  of  an  individual  subsystem  at  an 
individual  frequency  by  E0,  the  total  computational  effort  in  frequency  computations  for 
the  two  techniques  is 

E„  =  Xm£0  =  (MT0)-[(/2  (34) 

n=l 

for  the  direct  approach,  and 


11 


(35) 


££=S'».E.=S£or(/2-/,)-7',l, 

«  =  1  rt  =  l 

for  the  efficient  technique.  The  number  of  time  points  for  the  input  (nin  as  in  equation  (17)) 
is  the  same  for  both  techniques.  The  number  of  time  points  for  the  output  for  the  direct 
approach  is  given  by  equation  (18),  and  for  the  efficient  technique  is  given  by  the  same 
expression  with  T0  replaced  by  a  summation  over  all  of  {Tn}.  Since  the  Fourier  convolution 
computations  are,  in  general,  not  as  costly  computationally  as  the  running  of  the  steady- 
state  model,  they  have  little  impact  on  overall  computational  effectiveness.  (The  relative 
costs  of  the  convolutions  used  in  the  two  techniques  will  be  discussed.) 

From  the  computational  requirements  shown  in  equations  (34)  and  (35),  it  is  clear 
that  the  efficient  technique  is  appropriately  named  in  most  cases.  In  fact,  the  difference 
in  computational  effort  between  the  two  techniques  can  be  written  as 


AE  —  Eg  Eg 


V  N  \ 

- 

2t.-nt0 

to-/.) 

V"=1  ; 

<0, 


(36) 


with  equality  holding  only  when  Tn  =  T0  for  all  n.  This  expression  shows  the  improve¬ 
ment  in  using  the  computionally  efficient  technique  over  the  direct  approach  for  all  cases 
of  practical  interest.  Furthermore,  the  a  priori  determination  of  the  subsystem  T„  s  (or  an 
upper  bound  on  them)  is  usually  much  more  easily  obtained  than  an  a  priori  estimate  on 
the  upper  bound  of  To.  Thus,  in  practice,  the  value  of  To  that  must  be  used  is  very  large 
compared  to  the  true  value,  which  can  only  be  obtained  from  the  resulting  output.  This 
observation  will  become  evident  in  the  examples  found  in  the  following  sections. 

There  is  an  added  cost  for  the  efficient  technique  in  that  Fourier  convolutions  of 
the  form  of  equation  (29)  must  be  performed  for  each  subsystem.  However,  the  Fourier 
evaluation  of  the  incident  pulse  p{t)  only  needs  to  be  performed  once,  since  all  subsystems 
use  the  same  frequency  spacings.  Using  a  straightforward  inverse  Fourier  transform 
technique  for  the  convolutions  requires 


mn  *(2 mn)2  =0(mn3) 


(37) 


complex  multiplication  operations  for  each  subsystem  (other  simple  operations,  such  as 
addition,  are  not  included  due  to  their  relative  efficiency).  If  an  inverse  fast  Fourier 
transform  technique4  is  employed,  the  number  of  complex  multiplications  is 


mn  •  0(mn  log2  mn )  =  0(mn2  log2  mn ) . 


(38) 
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The  expressions  in  equations  (37)  and  (38)  give  the  number  of  multiplications 
required  to  create  the  interpolates  for  Fourier  series  that  are  evaluated  at  each  time  point  to 
arrive  at  pulse  responses  in  the  time  domain.  The  direct  technique  has  similar  expressions, 
but  is  only  performed  once  for  the  composite  system,  not  separately  for  each  of  the  N 
subsystems.  If  Fd  and  FE  represent  the  computational  cost  of  performing  the  required 
Fourier  convolutions  for  the  two  techniques,  the  computational  difference  in  the 
convolution  evaluations  can  be  written  as 

A F  =  Fe-Fd 

=  0(N  •  mn  log2  mn)-0(m log2  m).  ^ 


This  number  is  close  to  zero  in  most  applications  and  can  therefore  be  disregarded 
in  many  cases.  However,  when  there  are  many  subsystems  and  the  number  of  frequencies 
used  in  the  efficient  technique  is  not  appreciably  lower  than  the  direct  technique,  then  the 
computational  cost  of  the  Fourier  convolutions  should  be  taken  into  consideration. 

The  standard  technique  for  transient  pulse  response  modeling  for  multiple 
subsystems  is  transfer  function  addition.  This  technique  is  summarized  by  equations  (9) 
and  (20).  The  efficient  technique  that  has  been  developed  is  shown  by  equations  (25),  (28), 
(29),  and  (32),  and  has  been  shown  to  require  far  fewer  runs  of  the  steady-state  models  to 
maintain  the  same  degree  of  accuracy.  The  computational  effort  analysis  that  was 
performed  assumes  the  majority  of  the  computational  effort  involved  in  this  process  is  in 
running  a  steady-state  model  at  a  number  of  frequencies  for  each  of  the  subsystems 
(transfer  function  evaluation).  The  main  improvement  in  the  efficient  technique  is  in  the 
reduction  in  the  number  of  frequencies  that  are  required.  Two  examples  from  acoustics 
(highlight-based  acoustic  scattering  and  ray-based  acoustic  propagation)  will  be  shown  to 
illustrate  the  performance  improvement  of  the  efficient  technique. 
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ACOUSTIC  SCATTERING  MODELING  EXAMPLE 


An  example  of  a  modeling  task  that  can  make  use  of  the  technique  developed  in 
this  report  is  the  modeling  of  high-frequency  acoustic  scattering  using  multiple  highlight 
target  models.  In  these  applications,  the  acoustic  scattering  from  a  structure  is  represented 
as  a  composition  of  acoustic  scattering  from  a  number  of  discrete  features  (or  highlights) 
that  comprise  the  target.  These  individual  highlight  responses  can  then  be  linearly  super¬ 
imposed  to  create  the  composite  echo  return  from  the  structure.  Such  techniques  are 
common  in  high-frequency  scattering  where  theories  such  as  the  Geometric  Theory  of 
Diffraction  are  applicable.  In  these  cases,  there  is  a  natural  division  of  the  structure  into  a 
finite  number  of  dominant  highlights  (or  subsystems)  whose  linear  superposition  of 
responses  well  approximates  the  total  response.  The  main  goal  of  this  example  was  to 
show  the  utility  of  using  an  efficient  multiple  response  convolution  algorithm  as  opposed 
to  a  direct  technique.  In  that  sense,  it  is  assumed  that  the  direct  solution  is  correct,  that  this 
example  is  sufficiently  representative  of  the  general  acoustic  scattering  modeling  problem, 
and  that  the  results  drawn  from  it  may  be  extended  to  other  cases. 

For  the  example  in  this  report,  a  sample  target  comprised  of  15  highlights  was 
created.  This  target  was  modeled  in  a  steady-state  manner  using  the  high-frequency  target 
acoustic  model  (HFTAM)  to  obtain  the  complex  monostatic  response  at  a  variety  of 
frequencies  for  a  fixed  geometry.  The  HFTAM  target  model  is  used  in  high-frequency 
acoustic  scattering  modeling  applications  to  predict  the  response  (echo  return)  from  a 
submarine  structure  that  is  actively  ensonified.  Previous  applications  of  HFTAM  used 
a  narrowband  approximation  and  modeled  the  pulse  response  as  the  center  frequency 
response  only.  This  application  allows  for  the  computation  of  echo  returns  from  broadband 
pulses.  Two  geometries  were  used  in  this  study,  one  at  bow  aspect  and  the  other  at  near¬ 
beam  aspect  (80  degrees  off-bow).  In  the  near-beam  aspect  case,  all  of  the  highlights  will 
return  echoes  at  close  to  the  same  time,  creating  a  large  opportunity  for  coherent  addition. 
In  the  bow  aspect  case,  most  of  the  highlights  return  echoes  at  very  different  times, 
creating  very  little  opportunity  for  coherent  addition.  The  study  results  will  show  that  the 
efficient  computational  technique  performs  well  in  both  situations. 

For  this  example,  the  pulse  that  was  chosen  is  a  linear  frequency  modulator  (LFM) 
pulse  from  22  kHz  to  25  kHz  that  was  4  msec  in  duration.  The  time-bandwidth  product 
(TBP)  of  such  a  pulse  is  12,  which  is  broadband.  The  TBP  is  important  in  determining  the 
need  to  represent  a  pulse  with  nonzero  bandwidth  (broadband)  or  by  its  center  frequency 
alone  (narrowband).  If  a  pulse  has  a  large  bandwidth  but  a  very  short  time  duration,  there 
are  not  enough  cycles  of  energy  at  the  different  frequencies  to  make  a  broadband  repre¬ 
sentation  necessary.  Most  experts  agree  that  any  pulse  with  a  TBP  greater  than  7  is  a 
broadband  pulse  and  should  be  modeled  as  such.  Thus,  the  pulse  in  this  study  is  definitely 
a  broadband  pulse  and  should  be  modeled  with  all  of  the  necessary  component  frequencies. 
Note  that  the  bandwidth  used  in  computing  the  TBP  is  not  the  same  as  the  frequency  range 
used  in  the  previous  section  to  determine  frequency  spacing,  which  will  always  be  larger. 
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For  the  efficient  technique,  the  time  extent  for  each  of  the  subsystem  responses 
was  chosen  to  be  Tn  =  8  msec  (twice  the  input  pulse  length),  and  the  20-dB  down-point 
frequencies  used  in  the  modeling  were  20  kHz  and  27  kHz.  Based  on  these  numbers,  the 
frequency  spacing  necessary  for  pulse  response  modeling  is  125  Hz;  there  are,  therefore, 

57  frequencies  at  which  the  subsystem  steady-state  responses  (transfer  functions)  must  be 
computed.  For  the  direct  technique,  the  time  length  To  could  only  be  chosen  a  priori  to  be 
equal  to  twice  the  acoustic  length  of  the  target  plus  the  pulse  length,  or  141  msec  (the 
nominal  target  length  of  this  generic  structure  is  100  m).  This  is  a  definite  drawback  of  the 
direct  technique,  where  time  extent  is  often  chosen  to  be  greater  than  necessary  when  the 
results  of  the  model  output  are  not  available  ahead  of  time.  In  this  case,  the  direct 
technique  required  a  frequency  spacing  of  7  Hz;  there  were,  therefore,  1000  frequencies  at 
which  the  subsystem  steady-state  responses  were  computed.  This  is  a  17.9:1  difference  in 
the  number  of  target  model  runs  that  must  be  computed  for  the  two  techniques.  Based  on 
these  results,  the  difference  in  the  number  of  complex  multiplications  used  in  the  Fourier 
convolution  for  the  two  techniques  can  be  shown  (see  equation  (39))  to  be 

A  F  =  Fe-Fd 

=  0(4878)-  0(9966)  *0.  (40) 

It  can  be  seen  from  this  assessment  that  the  computational  improvement  of  transfer 
function  generation  for  the  subsystems  (57  versus  1000  frequencies)  of  the  efficient 
technique  is  quite  dramatic  and  presents  a  major  savings  in  run-time.  The  computational 
loads  for  Fourier  convolution  for  the  two  different  techniques  are  of  the  same  order  of 
magnitude  and  are  considered  unimportant  in  a  comparative  analysis. 

The  results  of  modeling  the  near-beam  (80  degrees  off-bow)  scattered  response 
of  this  pulse  using  HFTAM  and  the  two  pulse  response  techniques  are  shown  in  figure  1. 
From  this  time  series  representation,  it  is  clear  that  the  two  pulses  are  nearly  indistinguish¬ 
able.  The  many  additional  frequencies  of  the  direct  technique  were  needed  so  that  the 
summation  of  Fourier  modes  provides  the  zero  solution  (seen  from  0  to  60  msec  and  80  to 
150  msec  in  figure  1).  Since  the  efficient  technique  only  computes  a  solution  where  it  is 
expected  to  be  nonzero,  there  are  much  fewer  Fourier  modes  required  (much  less 
cancellation  of  modes  to  create  a  zero  solution  for  a  length  of  time).  A  good  measure  of  the 
matching  of  the  two  techniques  is  to  compute  the  spectrum  of  each  pulse  response  over  a 
region  of  time.  An  accurate  match  in  both  time  and  frequency  domains  will  guarantee  that 
most  signal  processing  algorithms  cannot  distinguish  between  the  two  returns.  For  this 
case,  a  time  region  encompassing  a  majority  of  the  returned  energy  (60  to  90  msec)  was 
chosen.  Figure  2  shows  a  plot  of  these  spectra  as  well  as  the  spectrum  of  the  incident  pulse. 
The  spectra  of  the  two  pulse  returns  show  that  the  direct  and  efficient  techniques  provide 
nearly  identical  solutions  in  the  region  of  highest  energy  concentration.  The  efficient 
technique  has  a  stronger  trailing  tail  off  of  the  main  spectrum,  caused  by  the  sharp  cutoff  to 
zero  in  a  finite  time  interval,  whereas  the  direct  technique  will  always  have  an  exponential 
decay  to  zero.  Furthermore,  note  that  the  spectra  of  the  returns  have  a  different  shape  than 
that  of  the  incident  pulse.  This  is  an  effect  of  both  the  broadband  nature  of  the  individual 
subsystem  returns  as  well  as  the  coherent  addition  of  multiple  subsystems. 
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The  results  of  modeling  the  bow  aspect  scattered  response  of  the  same  incident 
pulse  using  HFTAM  and  the  two  pulse  response  techniques  is  shown  in  figure  3.  Once 
again,  the  time  response  of  the  result  using  the  two  techniques  appear  to  be  nearly 
identical.  As  in  the  case  of  the  near-beam  response,  the  best  comparison  of  the  two  results 
is  made  by  looking  at  the  frequency  domain  response.  Since  there  are  distinct  regions  of 
time  where  the  energy  of  the  return  is  concentrated,  multiple  spectra  will  be  computed  for 
each  of  these  regions  of  time. 


Figure  3.  Comparison  of  Direct  and  Efficient  Echoes  for 
Scattering  at  Bow  Aspect 


Figures  4  and  5  show  two  of  the  regions  of  time,  from  0  to  20  msec  and  from  40  to 
60  msec,  respectively.  The  first  of  these  regions  (shown  in  figure  4)  contains  the  response 
from  only  a  single  highlight,  and  the  match  between  techniques  is  nearly  identical.  How¬ 
ever,  it  should  be  pointed  out  that  even  in  this  single  highlight  case,  the  spectral  shape  of 
the  return  is  different  than  that  of  the  incident  pulse.  This  is  the  primary  feature  of  broad¬ 
band  pulse  response  models.  Also,  the  tail  of  the  spectrum  for  the  efficient  technique  is 
larger  than  that  of  the  direct  technique.  As  before,  this  is  a  consequence  of  the  finite  time 
duration  of  the  efficient  technique’s  response.  The  spectra  in  figured  contain  the  returns 
from  multiple  highlights  that  are  superimposed  over  each  other  in  time.  The  spectra  match 
very  closely  and  contain  a  combination  of  multiple  highlight  features  (superposition)  as 
well  as  broadband  characteristics  of  each  of  the  highlights.  The  region  of  time  beyond  the 
primary  energy  does  have  some  small  amplitude  highlight  responses,  and  their  spectra  are 
shown  in  figure  6.  It  can  be  seen  that  the  two  techniques  vary  slightly  in  their  responses; 
however,  this  is  because  the  direct  technique  has  some  of  the  exponential  time  decay  of  the 
larger  amplitude  returns  corrupting  these  very  small  amplitude  returns.  Since  this  is  only 
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an  effect  on  the  very  smallest  of  returns  (relatively  small  within  any  given  pulse  return),  it 
is  not  a  problem  for  the  overall  echo  response.  Thus,  the  efficient  technique  matches  the 
direct  technique  result  very  well  at  a  greatly  reduced  computational  cost. 


Figure  4.  Spectrum  of  the  Echoes  for  Scattering  at  Bow  Aspect 
with  a  Window  from  0  to  20  msec 
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Figure  5.  Spectrum  of  the  Echoes  for  Scattering  at  Bow  Aspect 
with  a  Window  from  40  to  60  msec 


Figure  6.  Spectrum  of  the  Echoes  for  Scattering  at  Bow  Aspect 
with  a  Window  from  60  to  150  msec 


ACOUSTIC  PROPAGATION  MODELING  EXAMPLE 

For  the  case  of  acoustic  propagation  modeling,  the  Naval  Undersea  Warfare 
Center’s  comprehensive  acoustic  system  simulation  (CASS)  model  was  used  in 
conjunction  with  the  Gaussian  ray  acoustic  bundling  (GRAB)  propagation  model5  along 
with  Navy  standard  bottom  and  surface  bounce  models.  This  combination  of  models  has 
been  proven  to  be  effective  in  modeling  acoustic  reverberation  data.  The  main  goal  of 
this  study,  as  that  of  the  scattering  example,  was  to  show  the  utility  of  using  an  efficient 
multiple  response  convolution  algorithm  as  opposed  to  a  direct  technique.  In  that  sense,  it 
was  assumed  that  the  direct  solution  is  correct,  and  that  this  example  is  sufficiently 
representative  of  the  general  acoustic  propagation  modeling  problem  that  results  drawn 
from  it  may  be  extended  to  other  cases. 

The  runs  that  were  made  all  used  a  range  of  2  kyd  and  the  same  incident  pulse  as 
used  in  the  acoustic  scattering  example  (a  4-msec  LFM  from  22  to  25  kHz).  The  environ¬ 
mental  model  performs  a  ray  trace6  for  a  sequence  of  test  rays,  spaced  every  0.25  degrees 
in  the  vertical  aperture.  The  CASS  model  (with  GRAB  propagation  sub-models) 
computes  a  Gauss-weighted  averaging  of  these  rays  into  a  set  of  eigenrays,  whose 
combined  pressure  levels  mimic  the  true  pressure.  In  the  context  of  this  report,  the 
division  of  the  acoustic  propagation  problem  into  eigenrays  is  the  subsystem  breakdown 
that  will  be  used  in  the  efficient  computational  technique. 
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For  the  efficient  technique,  the  time  extent  for  each  of  the  subsystem  responses 
was  chosen  to  be  Tn  =  8  msec  (twice  the  input  pulse  length),  and  the  20-dB  down-point 
frequencies  used  in  the  modeling  were  20  kHz  and  27  kHz.  Based  on  these  numbers,  the 
frequency  spacing  necessary  for  pulse  response  modeling  is  125  Hz;  there  are,  therefore, 

57  frequencies  at  which  the  subsystem  steady-state  responses  (transfer  functions)  must  be 
computed.  There  were  also  nine  eigenrays  arriving  at  the  range  of  interest,  and,  thus,  nine 
subsystems.  For  the  direct  technique,  the  time  length  To  could  only  be  chosen  after  running 
the  model  once  at  a  nominal  frequency,  the  center  frequency  of  23.5  kHz  in  this  case,  and 
comparing  the  relative  time  delays  of  the  eigenrays.  The  difference  between  the  maximum 
and  minimum  eigenray  delay  times  was  54  msec;  adding  double  the  pulse  length  to  this 
makes  To  =  62  msec.  The  requirement  to  make  a  single  frequency  run  of  the  environmental 
model  before  choosing  To  is  a  distinct  drawback  to  the  direct  technique.  Based  on  this 
analysis,  the  direct  technique  required  a  frequency  spacing  of  16  Hz,  and  there  were, 
therefore,  438  frequencies  at  which  the  subsystem  steady-state  responses  were  computed. 
This  is  a  7.7:1  difference  in  the  number  of  target  model  runs  that  must  be  computed  for  the 
two  techniques.  Based  on  these  results,  the  difference  in  the  number  of  complex  multipli¬ 
cations  used  in  the  Fourier  convolution  for  the  two  techniques  can  be  shown  (see  equation 
(39))  to  be 

A  F  =  FE-FD 

=  0(2992)  -  0(3843) «  0.  (41) 

It  can  be  seen  from  this  assessment  that  the  computational  improvement  of  transfer 
function  generation  for  the  subsystems  (57  versus  438  frequencies)  of  the  efficient 
technique  is  dramatic  (although  not  as  much  as  in  the  scattering  example)  and  presents  a 
major  savings  in  runtime.  The  computational  loads  for  Fourier  convolution  for  the  two 
different  techniques  are  of  the  same  order  of  magnitude  and  are  considered  unimportant 
in  a  comparative  analysis. 

The  results  of  modeling  the  propagation  of  the  broadband  pulse  with  the  two 
techniques  are  shown  in  figure  7.  From  this  time  series  representation,  it  is  clear  that  the 
two  pulses  are  nearly  indistinguishable.  The  many  additional  frequencies  of  the  direct 
technique  were  needed  so  that  the  summation  of  Fourier  modes  provides  the  zero  solution. 
Since  the  efficient  technique  only  computes  a  solution  where  it  is  expected  to  be  nonzero, 
there  are  much  fewer  Fourier  modes  required  (much  less  cancellation  of  modes  to  creates  a 
zero  solution  for  a  length  of  time).  To  compare  the  two  time  series,  the  spectrum  of  each 
pulse  response  over  a  region  of  time  is  again  computed.  An  accurate  match  in  both  time 
and  frequency  domains  will  guarantee  most  signal  processing  algorithms  cannot 
distinguish  between  the  two  returns.  For  this  case,  two  time  regions  were  examined.  The 
first  time  region,  from  5  to  13  msec,  had  multiple  eigenrays  arriving  at  time  delays  that 
caused  the  output  pulses  to  overlap.  (Its  spectra  is  shown  in  figure  8.)  Both  the  efficient 
and  direct  techniques  have  nearly  identical  spectra  near  the  peak,  and  their  spectra  are 
different  from  the  incident  pulse,  thus,  showing  an  effect  of  the  superposition  of  multiple 
eigenrays.  The  second  time  region  examined  is  from  32  to  40  msec,  where  only  a  single 
eigenray  returned  with  a  time  delay  in  that  region.  These  spectra  are  shown  in  figure  9  and 
again  match  each  other  very  well.  In  this  case,  the  spectra  are  similar  in  shape  (although 
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not  magnitude)  to  the  incident  spectrum,  an  effect  of  the  single  eigenray.  For  both  of  these 
cases,  as  in  the  scattering  examples,  the  efficient  technique  has  a  stronger  trailing  tail  off  of 
the  main  spectrum,  caused  by  the  sharp  cutoff  to  zero  in  a  finite  time  interval,  whereas  the 
direct  technique  will  always  have  an  exponential  decay  to  zero. 


Figure  7.  Comparison  of  Direct  and  Efficient  Echoes 
for  Propagation  at  a  Range  of  2  kyd 


Figure  8.  Spectrum  of  the  Pulses  for  2-kyd  Propagation 
with  a  Window  from  5  to  13  msec 


Figure  9.  Spectrum  of  the  Pulses  for  2-kyd  Propagation  with  a 
Window  from  32  To  40  msec 


CONCLUSIONS 


This  report  has  presented  two  techniques  for  convolving  a  time  domain  pulse  with 
model-generated  transfer  functions  of  systems  composed  of  multiple  subsystems.  The 
first  technique  was  a  direct  method  whereby  the  individual  subsystem  transfer  functions 
were  summed  and  the  resulting  total  system  transfer  function  was  convolved  with  the 
pulse.  An  alternative  technique  of  convolving  individual  subsystems  with  the  pulse  and 
adding  the  results  together  in  the  (discretized)  time  domain  was  also  presented.  Details  of 
the  computational  implementation  of  the  technique,  including  frequency  step,  frequency 
range,  time  step,  etc.,  were  derived  from  standard  modeling  and  processing  principles. 

An  analysis  of  the  computational  requirements  of  the  two  techniques  showed  that 
the  subsystem  technique  provided  a  more  efficient  calculation  almost  all  of  the  time;  thus, 
it  was  referred  to  as  the  “efficient”  technique.  The  two  main  components  of  any  of  these 
model-based  convolution  techniques  are  running  the  model  a  number  of  times  to  generate 
the  transfer  function,  and  performing  the  numerical  Fourier  convolution.  It  was  shown  that 
the  efficient  technique  always  requires  less  effort  for  the  former,  and  the  effort  of  the  latter 
is  almost  always  lower  for  the  efficient  technique  (the  exception  is  when  there  are  many 
subsystems  that  all  contain  approximately  the  same  time  delay).  Furthermore,  the  efficient 
technique  was  shown  to  require  less  a  priori  information  than  the  direct  technique,  and, 
thus,  is  more  suited  to  a  completely  automated  modeling  framework. 
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The  two  techniques  were  each  applied  to  two  examples  of  Navy  interest:  acoustic 
scattering  from  a  submerged  structure  and  acoustic  propagation  in  the  ocean  medium. 
Both  cases  illustrated  that  the  efficient  technique  can  significantly  improve  computational 
speed  while  maintaining  accuracy.  Analysis  of  both  time  and  frequency  domain 
descriptions  of  the  output  from  the  two  techniques  showed  the  results  to  be  virtually 
indistinguishable  between  the  two  techniques,  a  conclusion  that  held  for  both  areas  of 
application. 

The  results  of  this  report  identify  a  computationally  efficient  technique  for 
convolution  of  signals  with  linear  systems  composed  of  multiple  linear  subsystems.  Two 
examples  of  Navy  application  showed  the  effectiveness  and  utility  of  the  technique.  It  is 
recommended  that  this  efficient  technique  be  applied  to  model-based  convolution  when¬ 
ever  there  are  multiple  subsystems  being  modeled.  The  technique  will  be  especially 
important  in  broadband  pulse  convolution  with  models  in  a  simulation  setting. 
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